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Abstract: The neutron spectrum unfolding by Bonner Spheres Spectrometer (BSS) is considered a 
complex multidimensional model, which requires complex mathematical methods to solve the first 
kind of Fredholm integral equation. In order to solve the problem of the Maximum Likelihood 
Expectation Maximization (MLEM) algorithm which is easy to suffer the pitfalls of local optima and 
the particle swarm optimization (PSO) algorithm which is easy to get unreasonable flight direction 
and step length of particles, which leads to the invalid iteration and affect efficiency and accuracy, an 
improved PSO-MLEM algorithm, combined of PSO and MLEM algorithm, is proposed for neutron 
spectrum unfolding. The dynamic acceleration factor is used to balance the ability of global and local 
search, and improves the convergence speed and accuracy of the algorithm. Firstly, the Monte Carlo 
method was used to simulated the BSS to obtain the response function and count rates of BSS. In the 
simulation of count rate, four reference spectra from the IAEA Technical Report Series No. 403 were 
used as input parameters of the Monte Carlo method. The PSO-MLEM algorithm was used to unfold 
the neutron spectrum of the simulated data and was verified by the difference of the unfolded spectrum 
to the reference spectrum. Finally, the *Cf neutron source was measured by BSS, and the PSO- 
MLEM algorithm was used to unfold the experimental neutron spectrum. Compared with maximum 
entropy deconvolution (MAXED), PSO and MLEM algorithm, the PSO-MLEM algorithm has fewer 
parameters, and automatically adjusts the dynamic acceleration factor to solve the problem of local 
optima. The convergence speed of the PSO-MLEM algorithm is 1.4 times and 3.1 times that of the 
MLEM and PSO algorithms. Compared with PSO, MLEM and MAXED, the correlation coefficients 
of PSO-MLEM algorithm increase by 33.1%, 33.5% and 1.9%, and the relative mean errors decreased 
by 98.2%, 97.8% and 67.4%. 
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1. Introduction 


The precise measurement of the neutron spectrum distribution in the radiation field is a basic 
process in nuclear physics research, radiation protection and monitoring, reactor safety and control, 
particle accelerator and other nuclear technology fields. During the process of neutron spectrum 
measurement, it is necessary to unfold the neutron spectrum. Because the neutron spectrum unfolding 
procedure is an underdetermined problem with numerous viable solutions, a specialized optimization 
algorithm must be used to find the optimal solution |!“ 

The neutron spectrum unfolding is a process of solving underdetermined equations, and at the 
earliest the least square method is used to unfold the neutron spectrum. However, the unfolded 
spectrum obtained through multiple iterations is quite different from the reference spectrum. However, 
by adding the prior information, the results could be improved to a certain extent ©!. Engl et al. used a 
regularization method to unfold the neutron spectrum, which constrained the solution by adding the 
prior information to improve the stability and continuity of the unfolded spectrum '!. Reginatto et al. 
used the NE213 scintillation detector to measure the neutron spectrum. Moreover, the authors 
combined the Bayesian analysis method with the maximum entropy principle to unfold the neutron 
spectrum and to achieving relatively high accuracy '”*!, Artificial intelligence algorithms were 
frequently utilized in the neutron spectrum unfolding Pl. Freeman et al. based on BSS used a genetic 
algorithm to unfold the neutron spectrum. This method used a computer-simulated "chromosome" to 
encode the neutron spectrum, and then the individuals with the best fitness was chosen as the final 
result. Compared with other neutron spectrum unfolding methods, the extracted results show that the 
genetic algorithm can obtain more accurate unfolded spectrum in the absence of prior information l°. 
Subsequently, Sharghi Ido et al. developed two different artificial neural network methods to unfold 
the neutron spectrum (multi-layer and single-layer). With this method, the human nervous system was 
simulated and the neural network model was trained with known inputs and outputs to adjust its 
internal parameters to appropriate values. Finally, Am-Be and ***Cf neutron sources were measured, 
and the unfolded spectrum have a good agreement !!"!_ M. Shahmohammadi Beni et al. proposed the 
Maximum Likelihood-Expectation Maximization algorithm to unfold the neutron spectrum, and 
unfolded the Am—Be neutron spectrum in the absence of prior information and got good results!'”!. 
The MLEM algorithm has strong convergence ability and high stability in the early stage of iteration, 
but the later iteration is easy to suffer the pitfalls of local optima. 

The PSO algorithm was developed by Eberhart and Kennedy in 1995 and is based on bird foraging 
behavior''*!, More specifically, the PSO begins with a randomly generated initial population (known 
as a particle swarm) that is randomly distributed as a candidate solution in the search space. The PSO 
contains the velocity and position for each particle (solution vector). The speed of each particle must 
be dynamically modified based on its own flight experience, as well as the flight experiences of the 
surrounding particles. On top of that, the position of each particle must be updated by using the cost 
function l“, PSO algorithm has the advantage of global convergence, strong randomness of search, 
and adjustable local search ability and global search ability in the whole iteration process. H. 
Shahabinejad et al. developed a novel neutron spectrum unfolding code (SDPSO) based on particle 
swarm optimization algorithm, and performed neutron spectrum unfolding for Am—Be neutron 
sources |'*!, The disadvantage is that the highly random search process is easy to get unreasonable 
flight direction and step length of particles, which leads to the invalid iteration and affect efficiency 
and accuracy. In this regard, many scholars have studied traditional PSO from the perspective of 
parameter setting, convergence and combining with other algorithms to improve algorithms. Some 
authors use the inertia weights to balance the global search capability and the local search capability 
[16-18] Some authors propose combination algorithms, which are based on the global search advantage 
of PSO algorithm, such as the combination of PSO and genetic algorithm!!”!, PSO and k-clustering 


algorithm?”!, 

In this paper, an improved PSO-MLEM algorithm, combined of PSO algorithm and MLEM 
algorithm, is proposed for neutron spectrum unfolding. The PSO-MLEM algorithm has the advantages 
of PSO and MLEM, which can avoid the pitfalls of local optima through the adjustment of dynamic 
acceleration factors, and improve the convergence speed and accuracy. 


2. Method and Materials 
2.1 Theoretical model of Bonner Spheres Spectrometer 

The Bonner Spheres Spectrometer (BSS) is considered a typical neutron spectrum measurement 
instrument that offers a wide energy range, high sensitivity, isotropic response, consistent neutron flux, 
and ambient dose equivalent. The neutron is moderated by a polyethylene shell in the BSS. The center 
uses a thermal neutron detector, which can be one of the following three types: *He, '°B and °Li. The 
He neutron proportional detector with the highest thermal neutron cross section is selected as the 
detector. For unfolding the neutron spectrum, the response function for each sphere of BSS is required, 
and is commonly calculated by the Monte Carlo method. The measurement principle can be expressed 
as the first kind of Fredholm integral equation P173: 


N; + & = fe" Ri(Ej)O(B AE i= 12,40; f= 12,-,m (1) 


where, N; is the count rate of the ith Bonner Sphere; e; is the measurement error of ith Bonner 
Sphere. Emax and Emin are the highest and lowest energy of the neutron; Ri(E;) is the Bonner 
Sphere response function corresponding to the jth energy interval; ®(£;) is the fluence rate of 
neutrons in the jth energy interval; n is the number of Bonner Spheres; m is the number of energy 
intervals. 

Since the R;(E;) and N; are discretely distributed, Eq. (1) can be expressed in the discrete form 
of Eq. (2): 


Ni +& = LR; t=1,2,--,n; jf =1,2,--,m (2) 
Eq. (2) can also be expressed in matrix form as follows: 
N+e=Ro (3) 


where, N = (N1, Na, ..., Nn)" > @ = (P1, P2 ---, Ø@m)", R is the response matrix of n x m. 

Usually, N and R are known quantities, and the number of Bonner spheres n is often less than the 
number of the energy interval m. The Eq. (3) belongs to the underdetermined equations. In this work, 
the PSO-MLEM algorithm was used to solve the underdetermined problem in neutron spectrum 
unfolding. 
2.2 Principle of PSO-MLEM algorithm 

This paper combines the PSO algorithm with the MLEM algorithm, and the dynamic acceleration 
factor is used to balance the global and local search ability, which improves the convergence speed 
and calculation accuracy of the PSO-MLEM algorithm. The MLEM algorithm is used to limit the 
search direction of particles during each iteration of the PSO-MLEM algorithm. Figure 1(a) depicts 
the particle optimization process in the PSO-MLEM algorithm. The MLEM algorithm guides the 
flight direction of the particle swarm to ensure that the particles always fly to the solution space of the 
MLEM algorithm. The theoretical position calculated by the MLEM algorithm replaces the position 
determined by the particle inertia, which effectively reduces the possibility of search stagnation and 
missing the optimal solution. At the same time, the displacement of the particles is determined by both 
the PSO and MLEM algorithms during each iteration of the PSO-MLEM algorithm. As shown in the 
dynamic acceleration factor transformation Eq. (8), the PSO-MLEM algorithm adopts a smaller cı 


and a larger c, in the early stage, so that the particles can quickly narrow the range of the solution 
space with the help the MLEM algorithm. With the number of iterations increases, MLEM algorithm 
is easy to suffer the pitfalls of local optima. Therefore, the PSO-MLEM algorithm gradually increases 
c, and decreases cz in the iterative process, which strengthens the local search ability of the 
algorithm and helps MLEM algorithm to jump out of the pitfalls of local optima. 

The following 6 steps further explain the general procedure of the PSO-MLEM method for 
neutron spectrum unfolding: 

1) Initialization of the particle swarm 

In the PSO-MLEM algorithm, the initial step is to define the number of particles, the number of 
particles is set to 600, and each particle is defined as a solution vector that must be optimized p = 
[P1 P2 +» Pm]. The solution vector g is the neutron spectrum to be optimized in the process of 
neutron spectrum unfolding, while the initial position and velocity of the particle are initialized by 
random numbers between 0 and 1. 

2) Calculation of the cost function 

The following cost function is specified to limit the search range of the PSO-MLEM algorithm: 


[IN cai—Nobsl|” lp-@reell” 
Cost = UH —— (4) 


[IN obsll llererll" 
where, Nops is the count rate of Bonner spheres; @ is the solution of each particle in each iteration; 
Pref 1s the reference spectrum that contains the a priori information; Na is calculated by using Eq. 
(5) as follows: 

Neal = Rg. (5) 

Through the reference spectrum Pref and the count rate Nops, the above cost function limits the 
search direction of the particle. The first term of the Cost function is the distance between Na and 
Nops» and the second term is the distance between the updated spectrum g and the reference spectrum 
Pref- When the minimum of Cost function is obtained, the @ is the result of the neutron spectrum 
unfolding. 

3) Update the local optima solution of a single particle (pbest) 

The solution of each particle is calculated in the current iteration process. Then the cost value of 
the current solution is compared with the cost value of the historical local optima solution(pbes?). If 
the current cost value is smaller than the historical cost value, the local optima solution is updated to 
the solution of the current particle. 

4) Update the global optima solution of all particles (gbest) 

The cost value of the local optima solution (pbest) is compared with the cost value of the global 
optima solution (gbest). If the cost value of the local optima solution is smaller than the cost value of 
the global optima solution, the global optima solution is updated to the local optima solution of the 
particle. 

5) Update the velocity and solution vector of the particles 

According to Eq. (6) and (7), the velocity v;; and solution gj; of the jth dimension of the ith 
particle are updated: 


vii(t+1)=7° (pbest;(t) = O) +r: (gbest(t) = O) (6) 
pij t +1) = c1: v(t +1) + c2: pp ™(t+1) (7) 


where, v; is the velocity vector of the ith particle, ø; is the solution vector of the ith particle, the 


pbest; is the historical local optima solution of the ith particle, and gbest is the global optima 
solution of all particles, r is a random number between [0, 1], t is the number of iterations, c is 
the acceleration factor between [0, 1], as shown in Eq. (8): 


Oh Tae (8) 


where, t is the current number of iterations, and Tmax is the maximum number of iterations. pif'*™ 


is calculated according to the MLEM algorithm using Eq. (9): 


n Ni 
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6) Judging convergence 

The search is terminated when the cost value of the global optima solution achieves the 
predefined iterative convergence condition or when the current number of iterations reaches the 
maximum number of iterations. If the above conditions are not met, steps 2)-6) are repeated. 

The algorithm flow chart summarizes the preceding 6 steps as follows Fig. 1(b): 
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Fig. 1. Graphical depiction of PSO-MLEM algorithm. (a) Schematic diagram of PSO-MLEM particle flight; (b) Flow 
chart of PSO-MLEM algorithm. 


2.3 Verification by Monte Carlo simulation 

In order to verify the accuracy and convergence of the PSO-MLEM algorithm, the Monte Carlo 
method was used to simulate the response function of the Bonner Spheres Spectrometer. And the 
count rates of Bonner spheres were simulated while four reference spectra from the IAEA Technical 
Report Series No. 403!*4! were used as the input spectra of the Monte Carlo method. The MXD_FC33 
program in the UMG unfolding software package V3.3, which was based on the method of maximum 
entropy deconvolution (MAXED), was used to unfold the neutron spectrum!**!, The reference 
spectrum, adding a random error, was set as the prior information of the PSO-MLEM algorithm and 
MAXED. Combined with the response function and the count rates of Bonner spheres obtained by 
simulation, the neutron spectrum was unfolded using PSO-MLEM, PSO, MLEM and MAXED, 
respectively. 

2.3.1 Bonner Spheres Spectrometer response function simulation 

The Bonner spheres model refers to the parameters used by the Physikalisch Technische 
Bundesanstalt (PTB). The Bonner spheres spectrometer consists of 15 polyethylene moderated 
spheres of various diameters, with diameters of 0, 2.5, 3, 3.5, 4, 4.5, 5, 6, 7, 8, 9.5, 10, 12, 15, 18 
inches. The response function has an energy range of 1x10-°~398 MeV, which is divided into 60 
energy groups °% Each sphere has a thermal neutron detector in the center. The center detector 
chooses the SP90 *He proportional counter manufactured by Centronics of the UK. The Bonner sphere 
model is shown in Fig.2(a): 

The response functions of Bonner spheres spectrometer refer to the counts caused by unit incident 
neutrons emitted from parallel source terms. The number of reactions between the neutrons entering 
the detector and *He must be recorded when the Monte Carlo method was used to simulate the 
response function. The response functions of Bonner spheres spectrometer to a single-energy parallel 


neutron beam with a diameter of d and incident energy of E is defined as follows?”!: 
M 
Ra(En) = gas (10) 


where, Ra(E) is the energy response (cm?) when the incident neutron energy is E, Mg represents 
the count rate of the Bonner sphere with diameter d, and ®(E) is the neutron fluence (cm) when 
the incident neutron energy is E. Figure 2(b) depicts the response functions of 15 Bonner spheres as 
simulated by the Monte Carlo method. 
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Fig. 2 (Color online) Bonner sphere model and Responses functions. (a) Bonner sphere model and *He proportional counter; 


(b) Responses functions of Bonner spheres spectrometer simulated by Monte Carlo method 


As Fig. 2(b), the Bonner sphere with a small diameter has a high response in the low energy range, 
and the Bonner sphere with a large diameter has a high response in the high energy range. In addition, 
as the diameter of the Bonner sphere increases, the maximum response gets narrower and moves 
towards higher energies. 

2.3.2 Validation of the unfolded spectrum 

As Fig.3(a) (a1), the Input-Spec1 is from the CERN-CEC high energy reference field facility, 40 
cm thick Fe shield in the IAEA Technical Report Series No. 403, the Input-Spec2 is from Booster- 
Synchrotron of the Stanford Synchrotron Radiation Laboratory (SSRL) SPEAR in the IAEA 
Technical Report Series No. 403, the Input-Spec3 is from CERN Pu-Be calibration spectra in the 
IAEA Technical Report Series No. 403 and the Input-Spec4 is from Final Focus Test Beam Facility 
(FFTB) in the IAEA Technical Report Series No. 403. In the Fig.3(a) (al), the Ref-Specl, Ref-Spec2, 
Ref-Spec3 and Ref-Spec4, which are the spectra from Input-Spec1 to Input-Spec2 adding the + 10% 
random error, are the prior information for the PSO-MLEM algorithm and MAXED. Figure 3(a) (a2) 
shows the count rates of Bonner spheres spectrometer obtained when the Input-Spec1 to Input-Spec4 
were used for the input spectrum of simulation. With the response function and the count rates of 
Bonner spheres spectrometer obtained by simulation, the PSO-MLEM, PSO, MLEM and MAXED 
were used to unfold the neutron spectrum. 

Figure 3(b1) - (e1) depict the results of the unfolded spectra by the PSO-MLEM, PSO, MLEM 
and MAXED with the prior information (Ref-Specl to Ref-Spec4). Figure 3(b2) - (e2) is the ratio 
relation between Naef/Nobs> Npso mLeMm/Nobs> Npso/Nobs> Nmrem/Nobs and Nvaxep/Nobs - 


Nops are the count rates of 15 Bonner spheres, Naef (Naef = RPinput) is the result of the response 


function times the input spectrum. Besides, Npso muem(Vpso mMLEmM = R@pso mrem) is the result of 
the response function times the unfolded spectrum of the PSO-MLEM algorithm, Npso(Npso = 
R@pso) is the result of the response function times the unfolded spectrum of the PSO algorithm, 


NmLem(NmLem = R@miem) is the result of the response function times the unfolded spectrum of the 
MLEM algorithm, and Nyaxgp(NMaxep = R@maxep) is the result of the response function times the 
unfolded spectrum of the MAXED. The closer the ratio is to 1, the smaller the difference between the 
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Fig. 3. Unfolded spectra of Monte Carlo simulation data. (a) Input spectrum and count rates of Bonner spheres (a1: The solid 
line is the input spectrum of Monte Carlo method, and the dashed line is the prior information with a random error; a2: The 
count rates of 15 Bonner spheres); (b) Unfolded spectrum of the Specl (b1: Comparison of PSO-MLEM, PSO, MLEM and 
MAXED unfolded spectra with the input spectrum; b2: The ratio of Npso mtem, Npso, Nmtem, Nmaxep and Naef to Nobs); 
(c) Unfolded spectrum of the Spec2 (c1: Comparison of PSO-MLEM, PSO, MLEM and MAXED unfolded spectra with the 
input spectrum; c2: The ratio of Npso mtem> Npso» Nutem,> Nmaxep and Naef to Nops); (d) Unfolded spectrum of the Spec3 
(d1: Comparison of PSO-MLEM, PSO, MLEM and MAXED unfolded spectra with the input spectrum; d2: The ratio of 
Npso mtem> Npso, Nutem> Nmaxen and Naef to Nops); (e) Unfolded spectrum of the Spec4 (e1: Comparison of PSO-MLEM, 
PSO, MLEM and MAXED unfolded spectra with the input spectrum; e2: The ratio of Npso mtem, Npsos Nutem> Nmaxep 
and Naef to Nobs) 

The Pearson correlation coefficient was utilized for the evaluation of the correlation between the 
unfolded neutron spectrum gy! and the input spectrum g'"P"t, and its value ranges from +1 and — 
1. The r can calculate as following equation: 


cal input cal input 


mpi pi -DP 2X; 


T pcalginput = F Ss (11) 
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The relative mean error (RME) was also used to calculate the difference of the unfolded spectrum 
from the input spectrum. The RME can calculate as following equation: 


1 input 
Pi 9; 
input 
Pi 


1 
RME gcai ginput = ei (12) 


, : f input 
where, m is the sample size, represents the number of energy intervals, and o” and pr P 


respectively represent the normalized value of the ith channel corresponding to the unfolded spectra 


and the input spectrum. For Eq. (11), when r=0, the two variables are uncorrelated, when r>0, there 
is a positive correlation, when r<0, there is a negative correlation, and the degree of correlation is as 
follows: very strong correlation (0.8—1.0), strong correlation (0.6—0.8), moderate correlation (0.4—0.6), 
weak correlation (0.2—0.4), very weak correlation or no correlation (0.0—0.2). For equation Eq. (12), 
under ideal conditions (gp! = gi"P¥t), RME is 0. 


Table1 Correlation coefficient (r) and relative mean error (RME) of PSO-MLEM, PSO, MLEM and MAXED unfolded 


spectra 
Unfolding Input spectrum 
rand RME oo a 
method Input-Spec1 Input-Spec2 Input-Spec3 Input-Spec4 
PSO- T peal input 1.0000 0.9989 0.9940 0.9991 
MLEM RME gcaiginput 0.0132 0.0203 0.0664 0.0195 
T peal pinput 0.9862 0.8123 0.7467 0.8815 
PSO 
RME gai ginput 1.6720 2.1514 3.7010 0.2062 
T pcalginput 0.9877 0.9481 0.7447 0.8655 
MLEM 
RME peal ginput 0.6627 1.0882 3.0150 0.1978 
T geal input 0.9917 0.9037 0.9753 0.9914 
MAXED 
RME gcal pinput 0.1305 0.3528 0.2035 0.0719 


As can be seen from Table 1 that even if the prior information with a random error from the input 
spectra of the Monte Carlo method is used in the PSO-MLEM algorithm, the correlation coefficients 
of the PSO-MLEM unfolded spectrum of the four neutron sources are all above the value of 0.99 and 
the relative mean error RME are all less than 0.07%, indicating that the unfolded spectrum by the 
PSO-MLEM has a strong correlation with the input spectrum of the Monte Carlo method. The 
correlation coefficient of the PSO and MLEM algorithm is significantly lower than the PSO-MLEM 
algorithm. The Input-Spec3 fluctuate more strongly, the results of PSO and MLEM algorithms are 
lower accuracy. The correlation coefficient of Spec3 is only 0.74 in PSO and MLEM algorithms, 
which is lower than the PSO-MLEM algorithm (0.99), and the relative mean error is also more than 
the PSO-MLEM algorithm. The result shows that the PSO-MLEM algorithm has better robustness 
when the neutron spectrum is complex. At the same time, compared with the unfolded spectrum by 
MAXED, the correlation coefficient of PSO-MLEM is increased by 1.9%, and the relative mean error 
is decreased by 67.4%. 

The time complexity of the algorithm is represented by Big O notation, and the time complexity 
of the PSO-MLEM algorithm is O(T*N), which is linearly related to the number of iterations T and 
population N of the algorithm. Fig. 4 shows a plot between the minimum cost of Eq. (5) and the 
number of iterations of the PSO-MLEM, PSO and MLEM algorithms, which indicate the convergence 
of algorithm. From the Fig.4, the minimum cost of the PSO-MLEM algorithm is less than that of the 
PSO and MLEM algorithms when the convergence condition is achieved. And when the convergence 
condition is achieved, the iterations number of the PSO-MLEM algorithm is 500, the iteration number 
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of the MLEM algorithm is 700, the iterations number of the PSO algorithm is 1550. The convergence 
speed of the PSO-MLEM algorithm is 1.4 times and 3.1 times higher than that of MLEM and PSO 
algorithms. 


—— PSO-MLEM 
=== PSO 
—— MLEM 


Cost 


0 
10° 10! 10? 10° 104 
Iterations 
Fig. 4 (Color online) A plot of the algorithm convergence between the minimum cost and the number of iterations of the PSO- 
MLEM, PSO and MLEM algorithms. 

3. Experiment 

The Bonner Spheres Spectrometer designed by the China Institute of Atomic Energy was used 
for experiment to further evaluate the viability of the algorithm '7*!, As Fig.5(a), the Bonner Spheres 
Spectrometer consists of one SP9 *He proportional counter and ten polyethylene spheres with 
diameters of 1, 2.5, 3, 3.5, 4, 5, 6, 8, 10, 12 inches. The *He neutron proportional counter as the central 
detector is wrapped in 10 polyethylene spheres. The neutron energy range of the response function is 
between 9.441*10"!° and 18.84 MeV, which is divided into 207 energy groups. The response function 
of BSS is shown in Fig.5(b). 
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(b) 
Fig. 5. Experimental Bonner sphere model and Responses functions. (a) Photograph of Bonner spheres spectrometer; (b) 


Responses functions of Bonner spheres spectrometer 


The BSS of the China Institute of Atomic Energy is used to measure the reference **Cf neutron 
source, and the count rates of BSS N = (N4, N3, ..., N11)” is obtained and used to unfold the neutron 
spectrum by the PSO-MLEM algorithm. The reference spectrum and unfolded spectrum by the PSO- 
MLEM algorithm are shown in in Fig. 6. 
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Fig. 6. "Cf neutron source unfolded spectrum (a: Comparison of PSO-MLEM unfolded spectrum with reference spectrum; 


b: Ratio of Neas Naef and Nops of each detector) 


Table2 Correlation coefficient (r) and relative mean error (RME) of unfolded spectrum 


Correlation coefficient (r) and relative mean error (RME) 
252Cf neutron 


r MRE peal pret 


palgref 
source 


0.9980 0.1718 


The Table 2 is the Correlation coefficient (r) and relative mean error (RME) of unfolded spectrum 
by PSO-MLEM algorithm to the reference spectrum. The error analysis of the unfolded spectrum 


shows that the correlation coefficient of the unfolded spectrum gy! is more than the 0.99. The 
relative mean error of the unfolded spectrum is 0.1718%. 
4. Conclusion 

In this work, to solve the problem of the MLEM algorithm which is easy to suffer the pitfalls of 
local optima and the problem of the PSO algorithm which is easy to get unreasonable flight direction 
and step length of particles, which leads to the invalid iteration and affect efficiency and accuracy, an 
improved PSO-MLEM algorithm was proposed for neutron spectrum unfolding. Based on the 
traditional PSO algorithm, PSO and MLEM algorithm are mixed, and the convergence speed and 
accuracy of the algorithm are improved by dynamic acceleration factor, and the algorithm will not 
suffer the pitfalls of the local optima. 

The Monte Carlo method was employed to simulate the response functions of the BBS and four 
reference spectra from the IAEA Technical Report Series No. 403 were used as the input spectra. The 
effects caused by prior information with random error and fluctuate more strongly of spectrum on the 
unfolded spectrum were investigated. The results show that the PSO-MLEM algorithm is not sensitive 
to the fluctuation of prior information, and the unfolded spectrum has a strong correlation with the 
input spectrum of Monte Carlo method. In addition, the PSO-MLEM algorithm has strong robustness 
to the neutron spectrum with large fluctuations. Besides, PSO-MLEM algorithm is used to unfold the 
experimental data of *°’Cf neutron source, and the result shows that the PSO-MLEM algorithm is 
competent for neutron spectrum unfolding. 


Author contributions 

All authors contributed to the study conception and design. Material preparation, data collection and analysis were 
performed by Hong-Fei Xiao, Qing-Xian Zhang, He-Yi Tan, Bin Shi, Jun Chen, Zhi-Qiang Cheng and Jian Zhang. 
The first draft of the manuscript was written by Hong-Fei Xiao, Qing-Xian Zhang and all authors commented on 
previous versions of the manuscript. All authors read and approved the final manuscript. 


Data Availability Statement 
The data that support the findings of this study are openly available in Science Data Bank at 


https://doi.org/10.57760/sciencedb.07741 and https://cstr.cn/31253.11.sciencedb.07741. 


References 


[1] G. Medkour Ishak-Boushaki, M. Allab, Comparison between standard unfolding and Bayesian 
methods in Bonner spheres neutron spectrometry. Nucl. Instrum. Meth. A. 689, 35-39 (2012). 
https://doi.org/10.1016/j.nima.2012.06.016 

[2] H. Song, Y. Yuan, T. Peng, et al., Optimization study on neutron spectrum unfolding based on 
the least-squares method. Nucl. Sci. Tech. 29, 118 (2018). https://doi.org/10.1007/s41365-018-0454- 
5 

[3] M. Shahmohammadi Beni, D. Krstic, D. Nikezic, Studies on unfolding energy spectra of 
neutrons using maximum -likelihood expectation—maximization method. Nucl. Sci. Tech. 30, 134 
(2019). https://doi.org/10.1007/s41365-019-0662-7 

[4] M. Reginatto, P. Goldhagen, MAXED, a computer code for maximum entropy deconvolution of 
multisphere neutron spectrometer data. Health Phys. 77, 579-583 


(1999). https://doi.org/10.1097/00004032-19991 1000-00012 


[5] F.G. Perey, Least-squares dosimetry unfolding: The program STAY'SL, (Oak Ridge National 
Laboratory/TM-6062, 1977). https://doi.org/10.2172/5443582 

[6] H.W. Engl, Discrepancy principles for Tikhonov regularization of ill-posed problems leading to 
optimal convergence rates. J. Optimiz. Theory App. 52, 209-215 (1987). 
https://doi.org/10.1007/BF00941281 

[7] M. Reginatto, P. Goldhagen, S. Neumann, Spectrum unfolding, sensitivity analysis and 
propagation of uncertainties with the maximum entropy deconvolution code MAXED. Nucl. 
Instrum. Meth. A. 476, 242-246 (2002). https://doi.org/10.1016/S0168-9002(01)01439-5 

[8] M. Reginatto, A. Zimbal, Bayesian and maximum entropy methods for fusion diagnostic 
measurements with compact neutron spectrometers. Rev. Sci. Instrum. 79, 23505 (2008). 
https://doi.org/10.1063/1.2841695 

[9] J. Wang, Y. Zhou, Z. Guo, et al., Neutron spectrum unfolding using three artificial intelligence 
optimization methods. Appl. Radiat. Isot. 147, 136-143 (2019). 
https://doi.org/10.1016/j.apradiso.2019.03.009 

[10] D.W. Freeman, D.R. Edwards, A.E. Bolon, Genetic algorithms - a new technique for solving 
the neutron spectrum unfolding problem. Nucl. Instrum. Meth. A. 425, 549-576 

(1999). https://doi.org/10.1016/S0168-9002(98)01427-2 

[11] A. Sharghi Ido, M.R. Bonyadi, G.R. Etaati, et al., Unfolding the neutron spectrum of a NE213 
scintillator using artificial neural networks. Appl. Radiat. Isotopes. 67, 1912-1918 

(2009). https://doi.org/10.1016/j.apradiso.2009.05.020 

[12] M. Shahmohammadi Beni, D. Krstic, D. Nikezic, et al., Studies on unfolding energy spectra of 
neutrons using maximum-likelihood expectation—maximization method. Nucl. Sci. Tech. 30, 134 
(2019). https://doi.org/10.1007/s41365-019-0662-7 

[13] R. Eberhart, J. Kennedy, A new optimizer using particle swarm theory. Proceedings of the sixth 
international symposium on micro machine and human science. IEEE 39-43 (1995). 

[14] H. Wang, W. Wang, Z. Wu, Particle swarm optimization with adaptive mutation for 
multimodal optimization. Appl. Math. Comput. 221, 296-305 (2013). 
https://doi.org/10.1016/j.amc.2013.06.074 

[15] H. Shahabinejad, M. Sohrabpour, A novel neutron energy spectrum unfolding code using 
particle swarm optimization. Radiat. Phys. Chem. 136, 9-16 (2017). 
https://doi.org/10.1016/j.radphyschem.2017.03.033 

[16] Y. Shi, R.C. Eberhart, Parameter selection in particle swarm optimization. International 
conference on evolutionary programming, Springer, 591—600 (1998). 

[17] L. Zhang, Y. Tang, C. Hua, et al., A new particle swarm optimization algorithm with adaptive 
inertia weight based on Bayesian techniques. Appl. Soft Comput. 28, 138-149 (2015). 
https://doi.org/10.1016/j.asoc.2014.11.018 

[18] I.K. Gupta, A. Choubey, S. Choubey, Particle swarm optimization with selective multiple 
inertia weights. 2017 8th International Conference on Computing, Communication and Networking 
Technologies (ICCCNT). IEEE Computer Society pp. 1—6 (2017). 

https: //doi. org/10.1109/iccent.2017.8204132 

[19] M.A. Memon, M.D. Siddique, S. Mekhilef et al., Asynchronous Particle Swarm Optimization- 
Genetic Algorithm (APSO-GA) Based Selective Harmonic Elimination in a Cascaded H-Bridge 
Multilevel Inverter. IEEE T. Ind. Electron. 69, 1477-1487 (2022). 
https://doi.org/10.1109/TIE.2021.3060645 


[20] J. Li, L. Zhang, X. Gu, A composite particle swarm optimization algorithm for hospital 
equipment management risk control optimization and prediction. J. Env. Public Health 2022, 1-9 
(2022). https://doi.org/10.1155/2022/5268887 

[21] T. Belgya, L. Szentmikl6si, Monte-Carlo calculated detector response functions to unfold 
radiative neutron capture spectra. Nucl. Instrum. Meth. A. 991, 165018 (2021). 
https://doi.org/10.1016/j.nima.2021.165018 

[22] V. Suman, P.K. Sarkar, Neutron spectrum unfolding using genetic algorithm in a Monte Carlo 
simulation. Nucl. Instrum. Meth. A. 737, 76—86 (2014). https://doi.org/10.1016/j.nima.2013.11.012 
[23] S. Maeda, H. Tomita, J. Kawarabayashi et al., Fundamental study on neutron spectrum 
unfolding using maximum entropy and maximum likelihood method. Prog. Nucl. Sci. Technol. 1, 
233-236 (2011). https://doi.org/10.15669/pnst. 1.233 

[24] International Atomic Energy Agency, Compendium of neutron spectra and detector responses 
for radiation protection purposes. Technical Reports Series No.403 (2001). 

[25] R. Li, J.B Yang, X.G Tuo, et al., Unfolding neutron spectra from water-pumping-injection 
multilayered concentric sphere neutron spectrometer using self-adaptive differential evolution 
algorithm. Nucl. Sci. Tech. 32, 26 (2021). https://doi.org/10.1007/s41365-021-00864-4 

[26] B. Wiegel, A.V. Alevra, B. Siebert, Calculations of the response functions of Bonner spheres 
with a spherical *He proportional counter using a realistic detector model. SCAN-9503183, (1994). 
[27] X. Jiang, J. Cao, C. Jiang, et al., Geant4 simulation of multi-sphere spectrometer response 
function and the detection of 4!Am—Be neutron spectrum. Nucl. Sci. Tech. 28, 174 (2017). 
https://doi.org/10.1007/s41365-017-0328-2 

[28] J. Chen, C.J. LI, W. Li, et al., Neutron spectra measurement of IHNI-I BNCT beam with multi- 
sphere spectrometer. Atomic Energy Science and Technology 48, 2127-2132 (2014). 
https://doi.org/10.7538/yzk.2014.48.11.2127 (in Chinese) 


